Load the necessary libraries
library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for effects plots in ggplot
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(DHARMa) # for assessing dispersion etc
library(lme4) # for lmer
library(lmerTest) # for degrees of freedom in lmer
library(glmmTMB) # for glmmTMB
library(performance) # for diagnostic plots
library(see) # for diagnostic plots
In an honours thesis from (1992), Mullens was investigating the ways that cane toads ( Bufo marinus ) respond to conditions of hypoxia. Toads show two different kinds of breathing patterns, lung or buccal, requiring them to be treated separately in the experiment. Her aim was to expose toads to a range of O2 concentrations, and record their breathing patterns, including parameters such as the expired volume for individual breaths. It was desirable to have around 8 replicates to compare the responses of the two breathing types, and the complication is that animals are expensive, and different individuals are likely to have different O2 profiles (leading to possibly reduced power). There are two main design options for this experiment;
Mullens decided to use the second option so as to reduce the number of animals required (on financial and ethical grounds). By selecting this option, she did not have a set of independent measurements for each oxygen concentration, by repeated measurements on each animal across the 8 oxygen concentrations.
Toad
Format of mullens.csv data file
| BREATH | TOAD | O2LEVEL | FREQBUC | SFREQBUC |
|---|---|---|---|---|
| lung | a | 0 | 10.6 | 3.256 |
| lung | a | 5 | 18.8 | 4.336 |
| lung | a | 10 | 17.4 | 4.171 |
| lung | a | 15 | 16.6 | 4.074 |
| ... | ... | ... | ... | ... |
| BREATH | Categorical listing of the breathing type treatment (buccal = buccal breathing toads, lung = lung breathing toads). This is the between subjects (plots) effect and applies to the whole toads (since a single toad can only be one breathing type - either lung or buccal). Equivalent to Factor A (between plots effect) in a split-plot design |
| TOAD | These are the subjects (equivalent to the plots in a split-plot design: Factor B). The letters in this variable represent the labels given to each individual toad. |
| O2LEVEL | 0 through to 50 represent the the different oxygen concentrations (0% to 50%). The different oxygen concentrations are equivalent to the within plot effects in a split-plot (Factor C). |
| FREQBUC | The frequency of buccal breathing - the response variable |
| SFREQBUC | Square root transformed frequency of buccal breathing - the response variable |
mullens <- read_csv("../public/data/mullens.csv", trim_ws = TRUE)
mullens %>% glimpse()
## Rows: 168
## Columns: 5
## $ BREATH <chr> "lung", "lung", "lung", "lung", "lung", "lung", "lung", "lung…
## $ TOAD <chr> "a", "a", "a", "a", "a", "a", "a", "a", "b", "b", "b", "b", "…
## $ O2LEVEL <dbl> 0, 5, 10, 15, 20, 30, 40, 50, 0, 5, 10, 15, 20, 30, 40, 50, 0…
## $ FREQBUC <dbl> 10.6, 18.8, 17.4, 16.6, 9.4, 11.4, 2.8, 4.4, 21.6, 17.4, 22.4…
## $ SFREQBUC <dbl> 3.255764, 4.335897, 4.171331, 4.074310, 3.065942, 3.376389, 1…
Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of copper, distance and their interaction on the number of number of worms. Area of the place segment was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual plates.
In this example, the individual TOADs are the random blocks. Each toad can only be of one BREATH type (they are either predominantly buccal breathers or predominantly lung breathers). Hence, BREATH is a between TOAD (block) effect. The frequency of buccal breathing of each TOAD was measured under eight different oxygen levels and thus, these represent the within block effect, as will the interaction between breathing type and oxygen level.
The response in this case is a little tricky since it is a proportion, but without the full original measurements. For example a frequency of 50% would indicate that half of the breaths taken by the toad during the monitoring phase were buccal breaths. Ideally, it would be good to have the actual counts (both the number of buccal breaths and the total number of breaths). That way, we could model the data against a binomial distribution.
As it is, we only have the proportion (as a percentage). Although we could model this against a beta distribution, this could be complicated if there are proportions of either 0 or 1 (100).
To help guide us through this and the other typical model assumptions, lets start by graphing the frequency of buccal breathing against breathing type and oxygen level. However, before we do, we need to make sure that all categorical variables are declared as factors (including the random effect). If we intend to model against a beta distribution, we will need a version of the response that is represented on a scale between 0 and 1 (but not include 0 or 1).
mullens <- mullens %>%
mutate(
BREATH = factor(BREATH),
TOAD = factor(TOAD),
pBUC = FREQBUC / 100,
pzBUC = ifelse(pBUC == 0, 0.01, pBUC)
)
So starting with the raw response.
ggplot(mullens, aes(y = FREQBUC, x = factor(O2LEVEL), color = BREATH)) +
geom_boxplot()
Conclusions:
If we intend to use a beta distribution, we could repeat the above with the scaled response. Since the rescaling maintains the same ranking and relative spacing, the plot will look the same as above, just with a different y-axis. However, our interpretation will change. Under a beta distribution (with logit link), we expect that the variance and mean will be related. We expect that the distributions will become more varied and more symmetrical as the expected mean shifts towards 0.5. Distributions approaching 0 and 1 will be more asymmetrical and smaller. This indeed does seem to be the case here.
Now lets explore the oxygen trends (separately for each breathing type).
ggplot(mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH)) +
geom_smooth() +
geom_point()
Conclusions:
ggplot(mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH)) +
geom_smooth() +
geom_point() +
facet_wrap(~ BREATH + TOAD, scales = "free")
## facet_grid(TOAD~BREATH)
Conclusions:
mullens.lmer1a <- lmer(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
data = mullens,
REML = TRUE
)
mullens.lmer1b <- lmer(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD),
data = mullens,
REML = TRUE
)
# OR
mullens.lmer1b <- update(mullens.lmer1a, . ~ . - (1 | TOAD) + (poly(O2LEVEL, 3) | TOAD))
mullens.allFit <- allFit(mullens.lmer1b)
## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]
mullens.allFit
## original model:
## SFREQBUC ~ BREATH + poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD) + BREATH:pol...
## data: mullens
## optimizers (7): bobyqa, Nelder_Mead, nlminbwrap, nmkbw, optimx.L-BFGS-B,nloptwrap.NLOPT_LN_N...
## differences in negative log-likelihoods:
## max= 0.987 ; std dev= 0.364
## Check which of the models are considered valid (OK)
is.OK <- sapply(mullens.allFit, is, "merMod")
is.OK
## bobyqa Nelder_Mead
## TRUE TRUE
## nlminbwrap nmkbw
## TRUE TRUE
## optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
## TRUE TRUE
## nloptwrap.NLOPT_LN_BOBYQA
## TRUE
AICc(mullens.lmer1a, mullens.lmer1b)
## mullens.lmer1c <- update(mullens.lmer1a, ~. -(1|TOAD) + (poly(O2LEVEL, 3)||TOAD))
## anova(mullens.lmer1c)
## mullens.lmer1d <- update(mullens.lmer1a, ~. + (0+poly(O2LEVEL, 3)||TOAD))
## anova(mullens.lmer1d)
Conclusions:
mullens.lmer1c <- lmer(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
data = mullens,
REML = FALSE
)
mullens.lmer1d <- update(mullens.lmer1c, . ~ . - BREATH:poly(O2LEVEL, 3))
AICc(mullens.lmer1c, mullens.lmer1d)
Conclusions:
mullens.glmmTMB1a <- glmmTMB(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
data = mullens,
family = gaussian(),
REML = TRUE
)
mullens.glmmTMB1b <- glmmTMB(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD),
data = mullens,
family = gaussian(),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
## OR
mullens.glmmTMB1b <- update(mullens.glmmTMB1a, ~ . - (1 | TOAD) + (poly(O2LEVEL, 3) | TOAD),
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
AICc(mullens.glmmTMB1a, mullens.glmmTMB1b)
Conclusions:
## mullens.glmmTMB1c <- glmmTMB(SFREQBUC ~ BREATH*poly(O2LEVEL, 3) + (1|TOAD), data=mullens,
## family=gaussian(), REML=FALSE)
## mullens.glmmTMB1d = update(mullens.glmmTMB1a, .~.-BREATH:poly(O2LEVEL, 3))
## AICc(mullens.glmmTMB1a, mullens.glmmTMB1b)
Conclusions:
mullens.glmmTMB2a <- glmmTMB(pzBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
data = mullens,
family = beta_family(link = "logit"),
REML = TRUE
)
mullens.glmmTMB2b <- glmmTMB(pzBUC ~ BREATH * poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD),
data = mullens,
family = beta_family(link = "logit"),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
## OR
mullens.glmmTMB2b <- update(mullens.glmmTMB2a, ~ . - (1 | TOAD) + (poly(O2LEVEL, 3) | TOAD),
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
AICc(mullens.glmmTMB2a, mullens.glmmTMB2b)
Conclusions:
## mullens.glmmTMB2a = glmmTMB(pzBUC ~ BREATH*poly(O2LEVEL, 3) + (1|TOAD), data=mullens,
## family=beta_family(link = "logit"), REML=FALSE)
## mullens.glmmTMB2b = update(mullens.glmmTMB2a, .~.-BREATH:poly(O2LEVEL, 3))
## AICc(mullens.glmmTMB2a, mullens.glmmTMB2b)
Conclusions:
plot_model(mullens.lmer1a, type = "diag")[-2] %>% plot_grid()
mullens.lmer1a %>% performance::check_model()
mullens.resid <- mullens.lmer1a %>% simulateResiduals(plot = TRUE)
plot_model(mullens.glmmTMB1b, type = "diag")[-2] %>% plot_grid()
mullens.glmmTMB1b %>% performance::check_model()
mullens.resid <- mullens.glmmTMB1b %>% simulateResiduals(plot = TRUE)
plot_model(mullens.glmmTMB2b, type = "diag")[-2] %>% plot_grid()
mullens.glmmTMB2b %>% performance::check_model()
## mullens.resid <- mullens.glmmTMB2a %>% simulateResiduals(plot=TRUE)
mullens.resid <- mullens.glmmTMB2b %>% simulateResiduals(plot = TRUE)
mullens.glmmTMB2b %>%
ggemmeans(~ O2LEVEL | BREATH) %>%
plot(add.data = TRUE)
mullens.lmer1a %>% summary()
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD)
## Data: mullens
##
## REML criterion at convergence: 453.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3421 -0.5127 0.0282 0.4435 3.4313
##
## Random effects:
## Groups Name Variance Std.Dev.
## TOAD (Intercept) 0.7733 0.8794
## Residual 0.7416 0.8612
## Number of obs: 168, groups: TOAD, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.77246 0.25810 19.00000 14.616 8.68e-12
## BREATHlung -1.00381 0.41817 19.00000 -2.400 0.0268
## poly(O2LEVEL, 3)1 -10.76355 1.09451 141.00000 -9.834 < 2e-16
## poly(O2LEVEL, 3)2 1.31085 1.09451 141.00000 1.198 0.2331
## poly(O2LEVEL, 3)3 -0.01429 1.09451 141.00000 -0.013 0.9896
## BREATHlung:poly(O2LEVEL, 3)1 13.03414 1.77330 141.00000 7.350 1.47e-11
## BREATHlung:poly(O2LEVEL, 3)2 -7.22946 1.77330 141.00000 -4.077 7.59e-05
## BREATHlung:poly(O2LEVEL, 3)3 2.75036 1.77330 141.00000 1.551 0.1231
##
## (Intercept) ***
## BREATHlung *
## poly(O2LEVEL, 3)1 ***
## poly(O2LEVEL, 3)2
## poly(O2LEVEL, 3)3
## BREATHlung:poly(O2LEVEL, 3)1 ***
## BREATHlung:poly(O2LEVEL, 3)2 ***
## BREATHlung:poly(O2LEVEL, 3)3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) BREATH p(O2LEVEL,3)1 p(O2LEVEL,3)2 p(O2LEVEL,3)3
## BREATHlung -0.617
## p(O2LEVEL,3)1 0.000 0.000
## p(O2LEVEL,3)2 0.000 0.000 0.000
## p(O2LEVEL,3)3 0.000 0.000 0.000 0.000
## BREATH:(O2LEVEL,3)1 0.000 0.000 -0.617 0.000 0.000
## BREATH:(O2LEVEL,3)2 0.000 0.000 0.000 -0.617 0.000
## BREATH:(O2LEVEL,3)3 0.000 0.000 0.000 0.000 -0.617
## BREATH:(O2LEVEL,3)1 BREATH:(O2LEVEL,3)2
## BREATHlung
## p(O2LEVEL,3)1
## p(O2LEVEL,3)2
## p(O2LEVEL,3)3
## BREATH:(O2LEVEL,3)1
## BREATH:(O2LEVEL,3)2 0.000
## BREATH:(O2LEVEL,3)3 0.000 0.000
mullens.lmer1a %>% tidy(conf.int = TRUE)
mullens.lmer1a %>%
tidy(conf.int = TRUE) %>%
kable()
| effect | group | term | estimate | std.error | statistic | df | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 3.7724598 | 0.258102 | 14.6161608 | 19 | 0.0000000 | 3.2322461 | 4.3126734 |
| fixed | NA | BREATHlung | -1.0038085 | 0.418173 | -2.4004623 | 19 | 0.0267793 | -1.8790545 | -0.1285624 |
| fixed | NA | poly(O2LEVEL, 3)1 | -10.7635498 | 1.094506 | -9.8341613 | 141 | 0.0000000 | -12.9273134 | -8.5997862 |
| fixed | NA | poly(O2LEVEL, 3)2 | 1.3108539 | 1.094506 | 1.1976670 | 141 | 0.2330562 | -0.8529097 | 3.4746175 |
| fixed | NA | poly(O2LEVEL, 3)3 | -0.0142947 | 1.094506 | -0.0130604 | 141 | 0.9895981 | -2.1780583 | 2.1494690 |
| fixed | NA | BREATHlung:poly(O2LEVEL, 3)1 | 13.0341442 | 1.773303 | 7.3502087 | 141 | 0.0000000 | 9.5284464 | 16.5398419 |
| fixed | NA | BREATHlung:poly(O2LEVEL, 3)2 | -7.2294601 | 1.773303 | -4.0768338 | 141 | 0.0000759 | -10.7351578 | -3.7237623 |
| fixed | NA | BREATHlung:poly(O2LEVEL, 3)3 | 2.7503615 | 1.773303 | 1.5509826 | 141 | 0.1231475 | -0.7553363 | 6.2560593 |
| ran_pars | TOAD | sd__(Intercept) | 0.8793850 | NA | NA | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 0.8611528 | NA | NA | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
mullens.lmer1a %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| SFREQBUC | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 3.77 | 0.26 | 3.26 – 4.28 | <0.001 |
| BREATH [lung] | -1.00 | 0.42 | -1.83 – -0.18 | 0.018 |
| O2LEVEL [1st degree] | -10.76 | 1.09 | -12.93 – -8.60 | <0.001 |
| O2LEVEL [2nd degree] | 1.31 | 1.09 | -0.85 – 3.47 | 0.233 |
| O2LEVEL [3rd degree] | -0.01 | 1.09 | -2.18 – 2.15 | 0.990 |
|
BREATH [lung] * O2LEVEL [1st degree] |
13.03 | 1.77 | 9.53 – 16.54 | <0.001 |
|
BREATH [lung] * O2LEVEL [2nd degree] |
-7.23 | 1.77 | -10.73 – -3.73 | <0.001 |
|
BREATH [lung] * O2LEVEL [3rd degree] |
2.75 | 1.77 | -0.75 – 6.25 | 0.123 |
| Random Effects | ||||
| σ2 | 0.74 | |||
| τ00 TOAD | 0.77 | |||
| ICC | 0.51 | |||
| N TOAD | 21 | |||
| Observations | 168 | |||
| Marginal R2 / Conditional R2 | 0.341 / 0.677 | |||
| AIC | 473.152 | |||
# with Satterthwaite degrees of freedom calculations
mullens.lmer1a %>% anova()
# OR with Keyward-Roger degrees of freedom calculations
mullens.lmer1a %>% anova(ddf = "Kenward-Roger")
mullens.glmmTMB1b %>% summary()
## Family: gaussian ( identity )
## Formula:
## SFREQBUC ~ BREATH + poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD) +
## BREATH:poly(O2LEVEL, 3)
## Data: mullens
##
## AIC BIC logLik deviance df.resid
## 452.6 511.9 -207.3 414.6 157
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## TOAD (Intercept) 0.8127 0.9015
## poly(O2LEVEL, 3)1 35.4538 5.9543 -0.11
## poly(O2LEVEL, 3)2 10.3349 3.2148 -0.30 -0.91
## poly(O2LEVEL, 3)3 2.3969 1.5482 0.13 -0.50 0.35
## Residual 0.4320 0.6573
## Number of obs: 168, groups: TOAD, 21
##
## Dispersion estimate for gaussian family (sigma^2): 0.432
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.77246 0.25821 14.610 < 2e-16 ***
## BREATHlung -1.00381 0.41835 -2.399 0.01642 *
## poly(O2LEVEL, 3)1 -10.76355 1.85070 -5.816 6.03e-09 ***
## poly(O2LEVEL, 3)2 1.31085 1.22182 1.073 0.28333
## poly(O2LEVEL, 3)3 -0.01429 0.93928 -0.015 0.98786
## BREATHlung:poly(O2LEVEL, 3)1 13.03414 2.99847 4.347 1.38e-05 ***
## BREATHlung:poly(O2LEVEL, 3)2 -7.22946 1.97958 -3.652 0.00026 ***
## BREATHlung:poly(O2LEVEL, 3)3 2.75036 1.52180 1.807 0.07071 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mullens.glmmTMB1b %>% tidy()
mullens.glmmTMB1b %>% tidy(effects = "fixed", conf.int = TRUE)
mullens.glmmTMB1b %>%
tidy(effects = "fixed", conf.int = TRUE) %>%
kable()
| effect | component | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|
| fixed | cond | (Intercept) | 3.7724598 | 0.2582114 | 14.6099634 | 0.0000000 | 3.2663746 | 4.2785449 |
| fixed | cond | BREATHlung | -1.0038085 | 0.4183504 | -2.3994445 | 0.0164200 | -1.8237601 | -0.1838568 |
| fixed | cond | poly(O2LEVEL, 3)1 | -10.7635498 | 1.8506966 | -5.8159452 | 0.0000000 | -14.3908484 | -7.1362512 |
| fixed | cond | poly(O2LEVEL, 3)2 | 1.3108539 | 1.2218219 | 1.0728682 | 0.2833303 | -1.0838731 | 3.7055808 |
| fixed | cond | poly(O2LEVEL, 3)3 | -0.0142947 | 0.9392759 | -0.0152188 | 0.9878576 | -1.8552416 | 1.8266523 |
| fixed | cond | BREATHlung:poly(O2LEVEL, 3)1 | 13.0341442 | 2.9984711 | 4.3469300 | 0.0000138 | 7.1572488 | 18.9110396 |
| fixed | cond | BREATHlung:poly(O2LEVEL, 3)2 | -7.2294601 | 1.9795777 | -3.6520213 | 0.0002602 | -11.1093612 | -3.3495590 |
| fixed | cond | BREATHlung:poly(O2LEVEL, 3)3 | 2.7503615 | 1.5218009 | 1.8073071 | 0.0707144 | -0.2323134 | 5.7330364 |
Conclusions:
# warning this is only appropriate for html output
mullens.glmmTMB1b %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| SFREQBUC | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 3.77 | 0.26 | 3.27 – 4.28 | <0.001 |
| BREATH [lung] | -1.00 | 0.42 | -1.82 – -0.18 | 0.016 |
| O2LEVEL [1st degree] | -10.76 | 1.85 | -14.39 – -7.14 | <0.001 |
| O2LEVEL [2nd degree] | 1.31 | 1.22 | -1.08 – 3.71 | 0.283 |
| O2LEVEL [3rd degree] | -0.01 | 0.94 | -1.86 – 1.83 | 0.988 |
|
BREATH [lung] * O2LEVEL [1st degree] |
13.03 | 3.00 | 7.16 – 18.91 | <0.001 |
|
BREATH [lung] * O2LEVEL [2nd degree] |
-7.23 | 1.98 | -11.11 – -3.35 | <0.001 |
|
BREATH [lung] * O2LEVEL [3rd degree] |
2.75 | 1.52 | -0.23 – 5.73 | 0.071 |
| Random Effects | ||||
| σ2 | 0.43 | |||
| τ00 TOAD | 0.81 | |||
| τ11 TOAD.poly(O2LEVEL, 3)1 | 35.45 | |||
| τ11 TOAD.poly(O2LEVEL, 3)2 | 10.33 | |||
| τ11 TOAD.poly(O2LEVEL, 3)3 | 2.40 | |||
| ρ01 | -0.11 | |||
| -0.30 | ||||
| 0.13 | ||||
| ICC | 0.72 | |||
| N TOAD | 21 | |||
| Observations | 168 | |||
| Marginal R2 / Conditional R2 | 0.338 / 0.813 | |||
| AIC | 452.576 | |||
mullens.glmmTMB2b %>% summary()
## Family: beta ( logit )
## Formula:
## pzBUC ~ BREATH + poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD) +
## BREATH:poly(O2LEVEL, 3)
## Data: mullens
##
## AIC BIC logLik deviance df.resid
## -493.0 -433.6 265.5 -531.0 157
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## TOAD (Intercept) 0.3140 0.5603
## poly(O2LEVEL, 3)1 11.9893 3.4626 0.11
## poly(O2LEVEL, 3)2 2.3761 1.5415 -0.26 -0.99
## poly(O2LEVEL, 3)3 0.8209 0.9060 -0.34 -0.85 0.90
## Number of obs: 168, groups: TOAD, 21
##
## Dispersion parameter for beta family (): 52.4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8176 0.1610 -11.291 < 2e-16 ***
## BREATHlung -0.5658 0.2622 -2.158 0.030912 *
## poly(O2LEVEL, 3)1 -6.7617 1.1184 -6.046 1.49e-09 ***
## poly(O2LEVEL, 3)2 0.5762 0.6893 0.836 0.403184
## poly(O2LEVEL, 3)3 0.1217 0.5932 0.205 0.837509
## BREATHlung:poly(O2LEVEL, 3)1 7.5618 1.8474 4.093 4.26e-05 ***
## BREATHlung:poly(O2LEVEL, 3)2 -4.0459 1.1650 -3.473 0.000515 ***
## BREATHlung:poly(O2LEVEL, 3)3 1.8213 1.0269 1.774 0.076133 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
mullens.glmmTMB2b %>% tidy(conf.int = TRUE)
mullens.glmmTMB2b %>% tidy(conf.int = TRUE, exponentiate = TRUE)
mullens.glmmTMB2b %>%
tidy(conf.int = TRUE, exponentiate = TRUE) %>%
kable()
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 0.1624188 | 0.0261464 | -11.2906065 | 0.0000000 | 0.1184701 | 2.226711e-01 |
| fixed | cond | NA | BREATHlung | 0.5678989 | 0.1488849 | -2.1582031 | 0.0309120 | 0.3397138 | 9.493555e-01 |
| fixed | cond | NA | poly(O2LEVEL, 3)1 | 0.0011573 | 0.0012944 | -6.0456042 | 0.0000000 | 0.0001292 | 1.036240e-02 |
| fixed | cond | NA | poly(O2LEVEL, 3)2 | 1.7793030 | 1.2264786 | 0.8359485 | 0.4031839 | 0.4608019 | 6.870456e+00 |
| fixed | cond | NA | poly(O2LEVEL, 3)3 | 1.1293592 | 0.6699145 | 0.2050814 | 0.8375085 | 0.3531144 | 3.612008e+00 |
| fixed | cond | NA | BREATHlung:poly(O2LEVEL, 3)1 | 1923.3855226 | 3553.3226762 | 4.0931655 | 0.0000426 | 51.4671290 | 7.187912e+04 |
| fixed | cond | NA | BREATHlung:poly(O2LEVEL, 3)2 | 0.0174948 | 0.0203821 | -3.4727187 | 0.0005152 | 0.0017833 | 1.716295e-01 |
| fixed | cond | NA | BREATHlung:poly(O2LEVEL, 3)3 | 6.1799931 | 6.3463423 | 1.7735771 | 0.0761331 | 0.8257990 | 4.624892e+01 |
| ran_pars | cond | TOAD | sd__(Intercept) | 0.5603453 | NA | NA | NA | 0.3974891 | 7.899258e-01 |
| ran_pars | cond | TOAD | sd__poly(O2LEVEL, 3)1 | 3.4625600 | NA | NA | NA | 2.2148406 | 5.413176e+00 |
| ran_pars | cond | TOAD | sd__poly(O2LEVEL, 3)2 | 1.5414544 | NA | NA | NA | 0.7802785 | 3.045171e+00 |
| ran_pars | cond | TOAD | sd__poly(O2LEVEL, 3)3 | 0.9060194 | NA | NA | NA | 0.1103516 | 7.438691e+00 |
| ran_pars | cond | TOAD | cor__(Intercept).poly(O2LEVEL, 3)1 | 0.1053372 | NA | NA | NA | -0.4110310 | 5.524276e-01 |
| ran_pars | cond | TOAD | cor__(Intercept).poly(O2LEVEL, 3)2 | -0.2594979 | NA | NA | NA | -0.2619112 | 2.622242e-01 |
| ran_pars | cond | TOAD | cor__(Intercept).poly(O2LEVEL, 3)3 | -0.3406095 | NA | NA | NA | -0.3386097 | 3.304468e-01 |
| ran_pars | cond | TOAD | cor__poly(O2LEVEL, 3)1.poly(O2LEVEL, 3)2 | -0.9862934 | NA | NA | NA | -0.7721392 | 9.492451e-01 |
| ran_pars | cond | TOAD | cor__poly(O2LEVEL, 3)1.poly(O2LEVEL, 3)3 | -0.8485818 | NA | NA | NA | -0.5798108 | 8.231511e-01 |
| ran_pars | cond | TOAD | cor__poly(O2LEVEL, 3)2.poly(O2LEVEL, 3)3 | 0.9000072 | NA | NA | NA | 0.8304543 | 8.480642e-01 |
# warning this is only appropriate for html output
mullens.glmmTMB2b %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| pzBUC | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 0.16 | 0.03 | 0.12 – 0.22 | <0.001 |
| BREATH [lung] | 0.57 | 0.15 | 0.34 – 0.95 | 0.031 |
| O2LEVEL [1st degree] | 0.00 | 0.00 | 0.00 – 0.01 | <0.001 |
| O2LEVEL [2nd degree] | 1.78 | 1.23 | 0.46 – 6.87 | 0.403 |
| O2LEVEL [3rd degree] | 1.13 | 0.67 | 0.35 – 3.61 | 0.838 |
|
BREATH [lung] * O2LEVEL [1st degree] |
1923.39 | 3553.32 | 51.47 – 71879.12 | <0.001 |
|
BREATH [lung] * O2LEVEL [2nd degree] |
0.02 | 0.02 | 0.00 – 0.17 | 0.001 |
|
BREATH [lung] * O2LEVEL [3rd degree] |
6.18 | 6.35 | 0.83 – 46.25 | 0.076 |
| Random Effects | ||||
| σ2 | 0.16 | |||
| τ00 TOAD | 0.31 | |||
| τ11 TOAD.poly(O2LEVEL, 3)1 | 11.99 | |||
| τ11 TOAD.poly(O2LEVEL, 3)2 | 2.38 | |||
| τ11 TOAD.poly(O2LEVEL, 3)3 | 0.82 | |||
| ρ01 | 0.11 | |||
| -0.26 | ||||
| -0.34 | ||||
| ICC | 0.71 | |||
| N TOAD | 21 | |||
| Observations | 168 | |||
| Marginal R2 / Conditional R2 | 0.334 / 0.809 | |||
| AIC | -492.981 | |||
The results so far suggest that the relationship between buccal breathing frequency and oxygen concentration is different for buccal breathers compared to lung breathers. We also know that there is strong evidence for a linear trend for the buccal breathers. It would therefore be interesting to specifically test whether there is evidence for either a linear or quadratic trend for the lung breathers.
mullens.lmer1a %>%
emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3) %>%
summary(infer = TRUE)
## OR
mullens.lmer1a %>% emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3, infer = c(TRUE, TRUE))
## degree = linear:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal -5.38e-02 1.33e-02 141 -8.01e-02 -2.75e-02 -4.048 0.0001
## lung -4.88e-03 1.69e-02 141 -3.84e-02 2.86e-02 -0.288 0.7740
##
## degree = quadratic:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal 4.45e-04 4.57e-04 141 -4.59e-04 1.35e-03 0.972 0.3326
## lung -2.67e-03 5.83e-04 141 -3.82e-03 -1.52e-03 -4.579 <.0001
##
## degree = cubic:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal -3.50e-07 2.69e-05 141 -5.36e-05 5.29e-05 -0.013 0.9896
## lung 6.73e-05 3.43e-05 141 -5.50e-07 1.35e-04 1.961 0.0518
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
newdata <- with(mullens, list(
O2LEVEL = modelr::seq_range(O2LEVEL, n = 100),
BREATH = levels(BREATH)
))
mullens.grid <- mullens.lmer1a %>%
emmeans(~ O2LEVEL | BREATH, at = newdata) %>%
as.data.frame()
mullens.grid %>%
group_by(BREATH) %>%
summarise(value = O2LEVEL[which.max(emmean)])
Conclusions:
mullens.glmmTMB1b %>%
emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3) %>%
summary(infer = TRUE)
## OR
emtrends(mullens.glmmTMB1b, specs = "BREATH", var = "O2LEVEL", max.degree = 3, infer = c(TRUE, TRUE))
## degree = linear:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal -5.38e-02 1.63e-02 157 -8.60e-02 -2.17e-02 -3.306 0.0012
## lung -4.88e-03 2.07e-02 157 -4.59e-02 3.61e-02 -0.235 0.8145
##
## degree = quadratic:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal 4.45e-04 4.48e-04 157 -4.40e-04 1.33e-03 0.993 0.3224
## lung -2.67e-03 5.71e-04 157 -3.80e-03 -1.54e-03 -4.675 <.0001
##
## degree = cubic:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal -3.50e-07 2.31e-05 157 -4.60e-05 4.53e-05 -0.015 0.9879
## lung 6.73e-05 2.95e-05 157 9.13e-06 1.26e-04 2.285 0.0236
##
## Confidence level used: 0.95
newdata <- with(mullens, list(
O2LEVEL = modelr::seq_range(O2LEVEL, n = 1000),
BREATH = levels(BREATH)
))
mullens.grid <- mullens.glmmTMB1b %>%
emmeans(~ O2LEVEL | BREATH, at = newdata) %>%
as.data.frame()
mullens.grid %>%
group_by(BREATH) %>%
summarise(value = O2LEVEL[which.max(emmean)])
Conclusions:
mullens.glmmTMB2b %>%
emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3) %>%
summary(infer = TRUE)
## OR
mullens.glmmTMB2b %>% emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3, infer = c(TRUE, TRUE))
## degree = linear:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal -3.46e-02 1.05e-02 157 -5.54e-02 -1.38e-02 -3.290 0.0012
## lung -9.07e-03 1.39e-02 157 -3.66e-02 1.85e-02 -0.651 0.5161
##
## degree = quadratic:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal 1.64e-04 2.13e-04 157 -2.57e-04 5.84e-04 0.769 0.4432
## lung -1.65e-03 3.52e-04 157 -2.34e-03 -9.53e-04 -4.683 <.0001
##
## degree = cubic:
## BREATH O2LEVEL.trend SE df lower.CL upper.CL t.ratio p.value
## buccal 2.99e-06 1.46e-05 157 -2.58e-05 3.18e-05 0.205 0.8378
## lung 4.78e-05 2.10e-05 157 6.32e-06 8.93e-05 2.276 0.0242
##
## Confidence level used: 0.95
## newdata <- with(mullens, list(O2LEVEL=seq(min(O2LEVEL), max(O2LEVEL), len=100),
## BREATH=levels(BREATH)))
mullens.grid <- with(mullens, list(
O2LEVEL = modelr::seq_range(O2LEVEL, n = 1000),
BREATH = levels(BREATH)
))
newdata <- emmeans(mullens.glmmTMB2b, ~ O2LEVEL | BREATH, at = mullens.grid) %>% as.data.frame()
newdata %>%
group_by(BREATH) %>%
summarise(value = O2LEVEL[which.max(emmean)])
## r.squaredGLMM(mullens.glmmTMB2a)
performance::r2_nakagawa(mullens.glmmTMB2b)
## # R2 for Mixed Models
##
## Conditional R2: 0.809
## Marginal R2: 0.334
Conclusions:
mullens.grid <- with(
mullens,
list(
BREATH = levels(BREATH),
O2LEVEL = modelr::seq_range(O2LEVEL, n = 100)
)
)
newdata <- mullens.lmer1a %>%
emmeans(~ O2LEVEL | BREATH, at = mullens.grid, type = "response") %>%
as.data.frame() %>%
mutate(across(c(emmean, lower.CL, upper.CL), function(x) x^2))
## mutate(emmean = emmean^2,
## lower.CL=lower.CL^2,
## upper.CL=upper.CL^2)
head(newdata)
ggplot() +
geom_ribbon(
data = newdata,
aes(
ymin = lower.CL, ymax = upper.CL,
x = O2LEVEL, fill = BREATH
), alpha = 0.3
) +
geom_line(
data = newdata,
aes(y = emmean, x = O2LEVEL, color = BREATH)
) +
theme_classic()
mullens.grid <- with(
mullens,
list(
BREATH = levels(BREATH),
O2LEVEL = modelr::seq_range(O2LEVEL, n = 100)
)
)
newdata <- mullens.glmmTMB1b %>%
emmeans(~ O2LEVEL | BREATH, at = mullens.grid, type = "response") %>%
as.data.frame() %>%
mutate(across(c(emmean, lower.CL, upper.CL), function(x) x^2))
## mutate(emmean = emmean^2,
## lower.CL=lower.CL^2,
## upper.CL=upper.CL^2)
head(newdata)
ggplot() +
geom_ribbon(
data = newdata,
aes(
ymin = lower.CL, ymax = upper.CL,
x = O2LEVEL, fill = BREATH
), alpha = 0.3
) +
geom_line(
data = newdata,
aes(y = emmean, x = O2LEVEL, color = BREATH)
) +
theme_classic()
mullens.grid <- with(
mullens,
list(
BREATH = levels(BREATH),
O2LEVEL = modelr::seq_range(O2LEVEL, n = 100)
)
)
newdata <- mullens.glmmTMB2b %>%
emmeans(~ O2LEVEL | BREATH, at = mullens.grid, type = "response") %>%
as.data.frame()
head(newdata)
ggplot() +
geom_ribbon(
data = newdata,
aes(
ymin = lower.CL, ymax = upper.CL,
x = O2LEVEL, fill = BREATH
), alpha = 0.3
) +
geom_line(
data = newdata,
aes(y = response, x = O2LEVEL, color = BREATH)
) +
scale_y_continuous("Buccal breathing rate", labels = function(x) 100 * x) +
theme_classic()
obs <- mullens %>%
mutate(
.fixed = predict(mullens.glmmTMB2b, re.form = NA),
.resid = residuals(mullens.glmmTMB2b),
PartialObs = plogis(.fixed + .resid)
)
ggplot() +
geom_point(data = mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH), alpha = 0.2) +
geom_point(data = obs, aes(y = PartialObs, x = O2LEVEL, color = BREATH)) +
geom_ribbon(
data = newdata,
aes(
ymin = lower.CL, ymax = upper.CL,
x = O2LEVEL, fill = BREATH
), alpha = 0.3
) +
geom_line(
data = newdata,
aes(y = response, x = O2LEVEL, color = BREATH)
) +
scale_y_continuous("Buccal breathing rate", labels = function(x) 100 * x) +
theme_classic()